/*
 * Copyright 2013 Sebastian Gesemann <s.gesemann@gmail.com>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * Simple implementation of the Fast Fourier Transform
 */

#include <cassert>
#include <cmath>
#include <algorithm>
#include "fft.hpp"

#define IS_POWER_OF_2(x) (((x)&((x)-1))==0)

namespace sg {

namespace {

using std::polar; // from <complex>
using std::swap;  // from <algorithm>

void shuffle(array_ref<complx> data)
{
	const auto n = data.size();
	assert(IS_POWER_OF_2(n));
	for (unsigned i=0, j=0; i<n; ++i) {
		if (i<j)
			swap(data[i],data[j]);

		// "reverse increment" j (with the bits' significance reversed)
		for (unsigned m=n/2; ((j^=m)&m)==0 && m>1; m/=2);
	}
}

void fft_core(array_ref<complx> data, int sign)
{
	const auto n = data.size();
	assert(IS_POWER_OF_2(n));
	complx cs = polar(real(1),sign*2*pi/n);
	for (unsigned x=n; x>=2;) {
		unsigned h = x/2;
		complx t (1,0);
		for (unsigned i=0; i<h; ++i) {
			for (unsigned l=0; l<n; l+=x) {
				complx& x = data[l + i];
				complx& y = data[l+h+i];
				complx  d = x-y;
				x = x + y;
				y = d * t;
			}
			t = t * cs;
		}
		x = h;
		cs = cs * cs;
	}
	shuffle(data);
}

} // anonymous namespace

void fft(array_ref<complx> data)
{
	fft_core(data,-1);
}

void ifft(array_ref<complx> data)
{
	const auto n = data.size();
	fft_core(data,1);
	const real scale = real(1)/n;
	for (unsigned i=0; i<n; ++i)
		data[i] *= scale;
}

} // namespace sg

